Python for Earth System Sciences

Remon Sadikni, Alexander Winkler

May 18th, 2018

---

Cartopy - an introduction

Geospatial data processing and projections

Overview

• Key features are the object-oriented projection definitions and its ability to transform points, lines, vectors, polygons and images between those projections
• Developed at the UK Met Office
• Basemap replacement (why still Basemap? $\Rightarrow$ Many people still use it; not all common projections yet in Cartopy)
• Publication quality geospatial projections

Source of this introduction: Cartopy Documentation

Getting Started

First, we have to import Matplotlib and Cartopy including its Coordinate Reference System.

## Import modules
# Matplotlib
%matplotlib notebook
import matplotlib.pyplot as plt
plt.ioff() # turn off interactive plotting mode

# Cartopy
import cartopy as cp

# Coordinate Reference System
import cartopy.crs as ccrs

Bad key "axes.color_cycle" on line 17 in
/Users/awi/.matplotlib/stylelib/mystyle.mplstyle.
You probably need to get an updated matplotlibrc file from
http://github.com/matplotlib/matplotlib/blob/master/matplotlibrc.template
or from the matplotlib source distribution

The simplest projection: Plate Carrée - equirectangular projection

We can access the individual projections from the imported cartopy coordinate reference system, e.g. ccrs.PlateCarree().

plt.figure(figsize=(10,5))

# Instantiate matplotlib axes with Cartopy Plate Carrée projection
ax = plt.axes(projection=ccrs.PlateCarree())

# Add coastlines
ax.coastlines();

plt.show()

Displaying available projections

Explore Cartopy's projections list or check

dir(ccrs)

# or hit TAB in an IPython Console

ccrs.<TAB>

directly in your console!

# Check for available projections
print dir(ccrs)

['ABCMeta', 'AlbersEqualArea', 'AzimuthalEquidistant', 'CRS', 'EuroPP', 'GOOGLE_MERCATOR', 'Geocentric', 'Geodetic', 'Geostationary', 'Globe', 'Gnomonic', 'InterruptedGoodeHomolosine', 'LambertAzimuthalEqualArea', 'LambertConformal', 'LambertCylindrical', 'Mercator', 'Miller', 'Mollweide', 'NearsidePerspective', 'NorthPolarStereo', 'OSGB', 'OSNI', 'Orthographic', 'PROJ4_VERSION', 'PlateCarree', 'Projection', 'Robinson', 'RotatedGeodetic', 'RotatedPole', 'Sinusoidal', 'SouthPolarStereo', 'Stereographic', 'TransverseMercator', 'UTM', 'WGS84_SEMIMAJOR_AXIS', 'WGS84_SEMIMINOR_AXIS', '_BoundaryPoint', '_CylindricalProjection', '_RectangularProjection', '_Satellite', '_WarpedRectangularProjection', '__builtins__', '__doc__', '__document_these__', '__file__', '__name__', '__package__', '__warningregistry__', '_ellipse_boundary', '_find_first_ge', 'absolute_import', 'abstractproperty', 'cartopy', 'division', 'epsg', 'math', 'np', 'prep', 'print_function', 'sgeom', 'six', 'warnings']

Let's choose an equal-area projection, e.g. Mollweide

# Mollweide projection with central longitude 180° W
plt.figure(figsize=(10,5))
ax = plt.axes(projection=ccrs.Mollweide(central_longitude=180))
ax.coastlines()

# Add stock background image
ax.stock_img()

plt.show()

Let's have a look at another one: Interrupted Goode homolosine projection

# Intantiate new figure
plt.figure(figsize=(10,5))
ax = plt.axes(projection=ccrs.InterruptedGoodeHomolosine()) 
ax.stock_img()

plt.show()

The feature interface

Cartopy feature interface relys on the comprehensive collection of free vector and raster map data provided by Natural Earth Data.

Natural Earth Data provides at 1:10m, 1:50m, and 1:110m scales.

---

Most commom features are pre-defined:

Feature	Description
cartopy.feature.BORDERS	Country boundaries
cartopy.feature.COASTLINE	Coastline, including major islands
cartopy.feature.LAKES	Natural and artificial lakes
cartopy.feature.LAND	Land polygons, including major islands
cartopy.feature.OCEAN	Ocean polygons
cartopy.feature.RIVERS	Single-line drainages, including lake centerlines

Example: Oceans and Lakes

plt.figure(figsize=(10,5))
ax = plt.axes(projection=ccrs.PlateCarree())

# Add oceans and high resolution shapefile of lakes
ax.add_feature(cp.feature.OCEAN)
ax.add_feature(cp.feature.LAKES.with_scale('10m'))

plt.show()

Your custom feature

---

Let's integrate our custom feature.

Example: I want to plot railroads.

1. Get download link of the data provided by Natural Earth Data: http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/cultural/ne_10m_railroads.zip

2. Figure out category (='cultural'), name (='railroads'), and scale (='10m').

3. Define your railroads feature (Cartopy gets the data for you):

   railroads = cp.feature.NaturalEarthFeature(category='cultural', 
                                               name='railroads', 
                                               scale='10m', 
                                               facecolor='none', # facecolor is set to 'none', because the polygons should not be color-filled
                                               color='blue') # color of the railroads
   

4. Add your railroads feature to your projection: ax.add_feature(railroads)

# Railroads feature example
plt.figure(figsize=(10,5))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.stock_img()

## Add a custom feature from Natural Earth Data: world's railroads
railroads = cp.feature.NaturalEarthFeature(category='cultural', 
                                              name='railroads', 
                                              scale='10m', 
                                              facecolor='none', # facecolor is set to 'none', because the polygons should not bet filled
                                              color='blue')
ax.add_feature(railroads)

<cartopy.mpl.feature_artist.FeatureArtist at 0x1809f56fd0>

plt.show()

Adding data to the map

---

Data are added to a projection just the same as they are added to a standard matplotlib figure.

Example: Plotting airlines between several cities

plt.figure(figsize=(10,5))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.stock_img()

# Hamburg coordinates
hh_lon, hh_lat = 10, 53

# Portland coordinates
pl_lon, pl_lat = -122, 45

# Sydney coordinates
sd_lon, sd_lat = 151, -33

# Plot connecting lines on a plane
plt.plot([pl_lon, hh_lon], [pl_lat, hh_lat],
         color='gray', linestyle='--', linewidth=3,)

plt.plot([sd_lon, hh_lon], [sd_lat, hh_lat],
         color='gray', linestyle='--', linewidth=3)

# Add labels for each cities
plt.text(pl_lon - 3, pl_lat - 8, 'Portland',
         horizontalalignment='right', fontsize=18)

plt.text(sd_lon - 3, sd_lat - 8, 'Sydney',
         horizontalalignment='right', fontsize=18)

plt.text(hh_lon + 3, hh_lat + 8, 'Hamburg',
         horizontalalignment='left', fontsize=18)

Text(13,61,'Hamburg')

plt.show()

Geocentric vs. Geodetic coordinate system

---

Geocentric

Cartersian coordinate system, where x, y, z are coordinates from the center of the Earth

Geodetic

Latitude/longitude coordinate system with spherical topology, geographical distance and coordinates in degrees

To get the actual airlines between two cities on Earth, the data has be transformed to spherical coordinates using the keyword transform and the method ccrs.Geodetic:

plt.figure(figsize=(10,5))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.stock_img()

# Plot connecting lines on a plane
plt.plot([pl_lon, hh_lon], [pl_lat, hh_lat],
         color='gray', linestyle='--', linewidth=3,)

plt.plot([sd_lon, hh_lon], [sd_lat, hh_lat],
         color='gray', linestyle='--', linewidth=3)

# Add labels for each cities
plt.text(pl_lon - 3, pl_lat - 8, 'Portland',
         horizontalalignment='right', fontsize=18)

plt.text(sd_lon - 3, sd_lat - 8, 'Sydney',
         horizontalalignment='right', fontsize=18)

plt.text(hh_lon + 3, hh_lat + 8, 'Hamburg',
         horizontalalignment='left', fontsize=18)

# Plot connecting lines on spheric coordinates
plt.plot([pl_lon, hh_lon], [pl_lat, hh_lat],
         color='blue', linewidth=3, marker='o',
         transform=ccrs.Geodetic())

plt.plot([sd_lon, hh_lon], [sd_lat, hh_lat],
         color='red', linewidth=3, marker='o',
         transform=ccrs.Geodetic())

[<matplotlib.lines.Line2D at 0x1809727e90>]

plt.show()

Understanding the keywords transform and projection

---

Cartopy's core concept is that the projection of your axes is independent of the coordinate system your data is defined in. The projection argument is used when creating plots and determines the projection of the resulting plot (i.e. what the plot looks like).

The transform argument to plotting functions tells Cartopy what coordinate system your data are defined in.

Cartopy always assumes that the your data are defined in the same coordinate system as the projection you use. So, we have to tell Cartopy the data's coordinate system $\Rightarrow$ most of the time lon/lat grids with cartesian coordinate system (ccrs.PlateCarree).

The safest thing to do is always provide the transform keyword regardless of the projection you are using, and avoid letting Cartopy make assumptions about your data’s coordinate system. Doing so allows you to choose any map projection for your plot and allow Cartopy to plot your data where it should be.

Adding gridded data

Read in netCDF data

# Load netCDF4 module
import netCDF4 as nc4

# NCEP reanalysis data: Monthly 2m air temperature
# On ICDC server:
# filename='/data/icdc/reanalyses/ncep_reanalysis1/DATA/2m_airtemp_monthly/air2m.mon.mean.nc'
filename='Data/air2m.mon.mean.nc'

# Open file and create file identifier
f = nc4.Dataset(filename,'r')

# Retrieve variables
lat = f.variables['lat'][:]
lon = f.variables['lon'][:]
air = f.variables['air'][:,:,:]

# Calculate climatologic mean for July
air_jul = air[6::12,:,:].mean(axis=0)

# Kelvin to Celsius
air_jul = air_jul - 273.15

Add the data on a projection

# Instantiate projection
plt.figure(figsize=(10,5))
ax = plt.axes(projection=ccrs.Mollweide(central_longitude=0))

# Plot the data in filled contours; Using transform=ccrs.PlateCarree() 
# we tell Cartopy the coordinate system the data is defined in
ax.contourf(lon, lat, air_jul[:,:], transform=ccrs.PlateCarree())

# Add coastlines
ax.coastlines()

# Set the extent of the Axes to the limits of the projection
ax.set_global()

plt.show()

Close the gap in cyclic coordinates

# Import add_cyclic_point methd
from cartopy.util import add_cyclic_point

# Close the gap in the data and the longitude array
air_jul, lon = add_cyclic_point(air_jul, coord=lon)

# Instantiate projection
plt.figure(figsize=(10,5))
ax = plt.axes(projection=ccrs.Mollweide(central_longitude=0))

# Plot the data in filled contours; Using transform=ccrs.PlateCarree() 
# we tell Cartopy the coordinate system the data is defined in
ax.contourf(lon, lat, air_jul[:,:], transform=ccrs.PlateCarree())

# Add coastlines
ax.coastlines()

# Set the extent of the Axes to the limits of the projection
ax.set_global()

plt.show()

Customize and polish the map

---

Example: 2m air temperature; climatologic mean for July; land without Antarctica

# Numpy
import numpy as np

# Import coordinates formating tool
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter

# Instantiate projection
plt.figure(figsize=(10,5))
ax = plt.axes(projection=ccrs.Miller(central_longitude=0))

# Mask ocean
ax.add_feature(cp.feature.OCEAN, zorder=1) # z order must be higher than the data

# Add coastlines
ax.coastlines()

# Add gridlines
ax.gridlines(color='k', linestyle='dashed')

# Cut off Antarctica
ax.set_extent([-180,180,-60,90]) # lon west, lon east, lat south, lat north

# Plot the data in filled contours, 17 levels from -36 to 36 deg C, specify colormap, extend colormap in both ways
p = ax.contourf(lon, lat, air_jul[:,:], 
                levels=np.linspace(-36, 36, 17), 
                extend='both', 
                transform=ccrs.PlateCarree(), 
                cmap=plt.get_cmap('RdBu_r'))

# Add colorbar
cbar = plt.colorbar(p, orientation='vertical', label='$^\circ$C', extend='both')

# Add longitude and latitude coordinates
ax.set_xticks([-180, -120, -60, 0, 60, 120, 180], crs=ccrs.PlateCarree())
ax.set_yticks([-60, -30, 0, 30, 60, 90], crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)

# Add title
plt.title('2m air temperature on land -- climatologic mean for July since 1948');

plt.show()

Another example: Quiver plot

---

Read in NCEP 10m wind data.

# Read ncep monthly wind data from ICDC
# on ICDC server:
# vwind = '/data/icdc/reanalyses/ncep_reanalysis1/DATA/2m_airtemp_monthly/vwnd10m.mon.mean.nc'
# uwind = '/data/icdc/reanalyses/ncep_reanalysis1/DATA/2m_airtemp_monthly/uwnd10m.mon.mean.nc'
fv = nc4.Dataset('Data/vwnd10m.mon.mean.nc','r')
fu = nc4.Dataset('Data/uwnd10m.mon.mean.nc','r')

# Retrieve variables
lat = fv.variables['lat'][:]
lon = fv.variables['lon'][:]
vwnd = fv.variables['vwnd'][:,:,:]
uwnd = fu.variables['uwnd'][:,:,:]

# Calc average winds for January and July
u_jan = uwnd[0::12,:,:].mean(axis=0)
u_jul = uwnd[6::12,:,:].mean(axis=0)

v_jan = vwnd[0::12,:,:].mean(axis=0)
v_jul = vwnd[6::12,:,:].mean(axis=0)

plt.figure(figsize=(10,5))
ax = plt.axes(projection=ccrs.Robinson())
ax.stock_img()

# Make streamplot
ax.quiver(lon, lat, u_jan, v_jan, transform=ccrs.PlateCarree(), regrid_shape=(len(lon)/4, len(lat)/4), width=0.001)

plt.title('Winds in January', fontsize=20)
plt.show()

plt.figure(figsize=(10,5))
ax = plt.axes(projection=ccrs.Robinson())
ax.stock_img()

# Make streamplot
ax.quiver(lon, lat, u_jul, v_jul, transform=ccrs.PlateCarree(), regrid_shape=(len(lon)/4, len(lat)/4), width=0.001)

plt.title('Winds in July', fontsize=20)
plt.show()

Of course, these plots need some polishing.

Plot data on curvilinear coordinates

import numpy as np
# Read in netcdf data on global tripolar grid: N deposition
ds = nc4.Dataset('TP04_ndepo_CMIP_NCAR_CCMI-1-0_gr_185001-185012-clim.nc','r')

# Retrieve variables
lat = ds.variables['lat'][:,:]
lon = ds.variables['lon'][:,:]
ndepo = ds.variables['ndepo'][:,:,:]

plt.figure(figsize=(10,3))

ax1 = plt.subplot(121, projection=ccrs.Robinson(central_longitude=0))
ax1.coastlines(resolution='110m', linewidth=0.5)

ax1.set_global()

## Plot using matplotlib contourf
p = ax1.contourf(lon, lat, ndepo[0,:,:],
               extend='max',
               levels=np.linspace(0.25e-11, 2e-11, 5),
               cmap = plt.get_cmap('Reds'),
               transform=ccrs.PlateCarree())

cbar = plt.colorbar(p, orientation='horizontal')
cbar.set_label('Ndepo')

ax1.set_title('Plot using contourf')

## Plot using matplotlib pcolormesh
ax2 = plt.subplot(122, projection=ccrs.Robinson(central_longitude=0))
ax2.coastlines(resolution='110m', linewidth=0.5)
ax2.set_global()

p = ax2.pcolormesh(lon, lat, ndepo[0,:,:],
               vmin=0.25e-11,
               vmax=2e-11,
               cmap = plt.get_cmap('Reds'),
               transform=ccrs.PlateCarree())

cbar = plt.colorbar(p, orientation='horizontal')
cbar.set_label('Ndepo')

ax2.set_title('Plot using pcolormesh')

Text(0.5,1,'Plot using pcolormesh')

plt.show()

Links

---

This was only a short introduction into Cartopy - there is a lot more to explore!

Start here:

• Outline $\Rightarrow$ Quickly find the reference documentation for known classes or functions of Cartopy.
• Gallery $\Rightarrow$ Many examples demonstrating some of the functionality of Cartopy, particularly its matplotlib interface.
• Projection List $\Rightarrow$ Browse the current list of projections.

Exercises:

---

1. Plot a global map of the climatologic mean of NCEP air temperature values in your favourite projection.
   • Add a colorbar.
   • Select a different colormap - here you find all available colormaps.
   • Add a grid.
2. Plot a map air temperature for the Mediterranean Sea (Longitudes = 8W -- 40E; Latitudes = 30N -- 46N).
   • Mask out the land.
   • Add a grid.
   • Add tick labels for longitude and latitudes.

3. Choose a feature on Natural Earth Data and visualize it in CartoPy (choose a projection you like).

4. Make a figure of the course of the river Elbe (or any river you like) - from its source to its outlet.